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Abstract 

We propose a simple mathematical model by applying Michaelis-Menton equations of enzyme kinetics 
to study the mutualistic interaction between the leaf cutter ant and its fungus garden at the early stage 
of colony expansion. We derive the sufficient conditions on the extinction and coexistence of these two 
species. In addition, we give a region of initial condition that leads to the extinction of two species when 
the model has an interior attractor. Our global analysis indicates that the division of labor by workers 
ants and initial conditions are two important factors that determine whether leaf cutter ants colonies and 
their fungus garden survive and grow can exist or not. We validate the model by doing the comparing 
between model simulations and data on fungal and ant colony growth rates under laboratory conditions. 
Wc perform sensitive analysis and parameter estimation of the model based on the experimental data 
to gain more biological insights on the ecological interactions between leaf cutter ants and their fungus 
garden. Finally, we give conclusions and discuss potential future work. 

Key words: Obligate Mutualism, Leaf-cutter ants. Garden Fungus, Functional/Numerical Response, 
Population Dynamics, Coexistence, Extinction, Multiple Attractors, Division of Labor, Parameter 
Estimation, Sensitive Analysis 



1. Introduction 

Mutualistic interactions, although ubiquitous in nature, are not well understood theoretically (Boucher 
1985; Herre and Bruna 1999; Hoeksema and Schwartz 2000; Holland and DeAngelis 2002; Neuhauser 
and Fargionc 2004). Mathematical modeling of mutualisms that correspond well with natural observed 
population dynamics have historically been difficult to formulate (Heithaus, Culver and Beattle 1980), 
especially for obligate mutualisms. While there are many obligate mutualisms found in nature, e.g., 
senita cacti and scnita moths, coral and zooxanthellae, plant and mycorrhizal fungi, termite and proto- 
zoa (Holland and DeAngelis 2010), there is little work done in mathematical modeling of this topic. Many 
two-species mutualism models come from modified Lotka-Volterra equations, which have been applied 
to a variety of ecological interactions, such as plant-pollinator interactions (Sobcron and Martinez del 
Rio 1981; Wells 1983; Holland and DeAngelis 2002) and legume rhizobium interactions (Vandermeer and 
Boucher 1978; Simms and Taylor 2002; West et al. 2002). However, these models predict unstable pop- 
ulation dynamics that do not match the dynamics actually observed in nature (Holland and DeAngelis 
2010). Recently, Holland and DeAngelis (2010) proposed a general framework using a consumer-resource 
approach to model the density-dependent population dynamics of mutualism. In this article, we apply 
Holland and DeAngelis' approach to derive a new and simple mathematical model of the population 
dynamics of leaf cutter ants and their fungus garden at the early colony stage, that can be validated by 
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the experimental data. This model is unique bcx;ause it ineludes the division of labor within the eolony 
and the particular colony stage; incorporation of these behavioral and life-history components make it 
distinct from other obligate mutualism models such as plant-pollinator interactions. 

This article focuses on a species of leaf-cutter ants; these ants are fungus farmers that harvest leaves 
and use them to cultivate their own food, a type of fungus, in underground gardens. The fungus that 
is grown by the adults feeds the ants' larvae, and the ants feed off of the fruit of fungus. Thus, leaf- 
cutter ants and their fimgus garden form an obligate mutualistic relationship, in which the increasing 
population of ants is due to consumption of the fungus while the increasing population of fungus is 
due to the agricultural services provided by the ants. The interaction between ants and fungus can be 
categorized as a consumer- resource mutualism according to Holland and DeAngelis's study (2010). 

Leaf-cutter ants of the genera Atta and Acromyrmex are among the most prevalent herbivores of the 
Neotropics, consuming far more vegetation than any other group of animals with comparable taxonomic 
diversity (Fowler et al. 1989; HoUdobler and Wilson 1990; Vieira-Neto and Vasconcelos 2010). In 
addition to their impact on plant communities, leaf-cutter ants can alter the spatio-temporal dynamics 
of carbon stocks, nutrient availability, susceptibility to fire, and other ecosystem properties through 
transferring tons of plant biomass below ground (Costa et al. 2008). Many leaf-cutter ant species increase 
in abundance following natural or anthropogenic disturbances in vegetation (reviewed in Wirth et al. 
2008). However, even though the harvesting behavior of ants is well-characterized, how this translates 
into the growth and production efficiency of the fungus and ant population is not well understood (Wirth 
2003). A mathematical treatment of the population dynamics between leaf-cutter ants and their fungus 
garden can be a first approach to explore these mechanisms. 

It is important to note that multiple mechanisms can potentially affect the performance of individuals 
and ultimately, regulate natural populations (Murdoch 1994), and that different mechanisms may operate 
at different life stages (Wilbur 1980; Jonsson and Ebenman 2001; Vieira-Neto and Vasconcelos 2010). 
For instance, colonies of most ant species are founded by solitary queens and must pass through a series 
of developmental stages before reaching maturity (reviewed in HoUdobler and Wilson 1990). Mortality 
is typically much higher for colonies in the early stages and therefore the forces acting on these colonies 
are likely to differ from those affecting older colonies (Wetterer 1994; Cole 2009). Studies examining 
variation across different stages of colony development are relatively rare and most have focused on 
ontogenetic changes in individual ant morphology and behavior (Tschinkel 1988; Wetterer 1994) rather 
than on factors affecting colony survival or growth. In this article we fociis on the ergonomic growth stage 
starting from when the first brood of workers reaches the adult stage. We study the population dynamics 
of leaf-cutter ants and their fungus garden during this critical stage by using a simple mathematical model 
to explore potential sources of colony failure and specific interactions that lead to different benefits and 
costs for colony growth. 

The model is unique also in that it incorporates behavioral effects, by considering the role of division 

of labor in colony growth, in particular the allocation of workers to different tasks that can positively or 
negatively affect fungal growth. Division of labor is one of the most basic and widely studied aspects of 
colony behavior in social insects, and addresses how individual worker behavior integrates into colony- 
level task organization (HoUdobler and Wilson 1990; Beshers and Fewell 2001). However, models of 
division of labor have not previously been integrated with effects on colony growth. In this paper, we 
adapted a model for incipient colony growth based on simple density-dependant ant growth and death 
rates coupled with a fungus growth model which can be described by a generalized Michaelis-Menton 
equations of enzyme kinetics model. The main purposes of this article are three-fold: 

1. Model functional response/numerical response based on ecological properties of leaf cutter ants 
and their fungus garden. 

2. Explore how the division of labor and initial conditions can be key factors that determine the 

successful colony expansion at its early stage. 

3. Validate the model with experimental data and perform parameter estimations and sensitivity 
analysis to understand the effects of parameters and initial condition on the model outcomes. 
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The rest of this article is organized as fohows: In section 2, we introduce the biological background 
of leaf-cutter ants and their fungus garden, and we formulate a simple mathematical model based on 
ecological assumptions that are supported by data and literature. In section 3, we perform mathematical 
analyses of the proposed model: We derive the sufficient conditions for the extinction and coexistence of 
the two species and give a region of initial conditions that leads to the extinction of ants and fungus with 
a model also containing an interior attractor (Theorem 3.1 and 3.2 1. These global analyses indicate that 
the division of labor by ants and initial conditions are two important factors in determining whether leaf 
cutter ants and their fungus garden can coexist or not (Theorem 3.3 and Corollary 3.1 ). In section 4, we 
compare simulations to data on growth rates for laboratory leaf cutter ant colonies, perform sensitivity 
analysis and parameter estimations for all the parameters and initial conditions around the nominal 
value. The study suggests not only that the fit of our model to data is significantly accurate, but also 
that the model can provide parameter values that are difficult to measure in the experiments. In the 
last section, we summarize our results and discuss future work. 



2. Biological Background and model formulation 

The life cycle of the leaf-cutter ant consists of four stages: egg, larva, pupa, and adult. Fertilized 
eggs produce female ants (queens and workers); unfertilized eggs produce male ants. The period from 
egg through adult stage usually lasts from 6 to 10 weeks. The average worker leaf-cutter ant lives from 
4-6 months in the lab (Howard, Henneman,Cronin, Fox and Hormig 1996), while the life span of workers 
in the field is shorter and varies by season. The queen is the largest of the ants and her life span can be 
10 years or more (Smith and Read 1997). 

Incipient colonies are started when mature colonies produce new winged females, queens, which will 
mate with males and then travel to a new location to start a new colony. They bring a small piece 
of fungus with them to establish a new garden in their new colony. After mating, queens construct 
nests by excavating in the soil. Within a day of nest initiation, the queen lays lays eggs, and uses her 
body's energy reserves to single-handedly raise her first brood, which consists entirely of workers. The 
first workers emerge 6-9 weeks after the eggs are laid (Holldobler and Wilson 1994). These workers 
enlarge the nest, feed the queen and larvae, tend to the fungus garden and hunt for leaves for the fungus 
substrate. Once the first workers emerge, the queens only task is to produce eggs. This transition begins 
the ergonomic growth phase of the colony, which may last for years until the colony matures. This 
growth stage is a critical time in the life of the colony, because the number of workers is low relative to 
the number of tasks and total work effort needed for colony growth (Brown, Bot and Hart 2006). 

Our study in this section aims to model the obligate mutualism interaction between leaf cutter ants 
and fungus at the incipient stage of the colony: From the time of first brood production to about 29 
weeks of age. During these early weeks, workers must allocate effort across a number of different tasks for 
colonies to survive and grow. The workers of established colonies, depending on maturity, may perform 
20-30 tasks (Wilson 1983). Incipient colonies perform considerably fewer; of these, feeding fungus to 
the larvae, tending the fungus garden, and collecting leaves for the fungus substrate require the bulk 
of worker activity budgets, which translate into energy expenditure of workers (Julian and Fewell 2004; 
Fewell unpublished data ). These tasks can be generalized into two categories; time spent outside the 
colony for collecting and processing leaves and time spent inside the colony for tending and cleaning 
the fungus garden and taking care of queens and larvae. For simplicity of modeling, we quantify time 
expenditure for these tasks as the biomass of ant. 

Let A(t) be the total biomass of ants including workers, larvae, pupae and eggs at time t, where pA 
( < p < 1) is the biomass of workers and (1 — p)A is the biomass of the remaining ants. Let F{t) be 
the total biomass of the fungus at time t. Ecological assumptions of the interaction between workers 
and fungus at the early stage of colony are as follows: 

Al: Assume that each worker has a fixed ratio of the energy spent outside the colony to the energy 
spent inside the colony which is This assumption is equivalent to a situation where workers 
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with a population of qpA (0 < q < 1) collect leaves, and the rest of workers (1 — q)pA tend the 
fungus garden and take care of queen ants as well as larvae. 

A2: The ants' population increases as the queen, larvae and adult ants feed on fungus. Thus we can 
assume that the numerical response function for ants is the HoUing Type I function, i.e., fungus 
biomass F multiplied by a constant number r^. In addition, we assume that ants suffer from 
density-dependent mortality due to energy consumed by foraging for leaves and taking care of 
the larvae and fungus garden, which will modify population growth through density-dependent 
self-limitation (Holland and DeAngelis 2010). Therefore, the population dynamics of ants can be 
described as follows: 

HA 

= {raF-d,A)A (1) 

where is a parameter that measures the maximum growth rate of ants and da is the mortality 
rate of ants. 

A3: The leaf-cutter ant mutualism is unique because the workers perform specific tasks to maintain the 
life of the fungus. The population of fungus can increase only if: there are qpA workers bringing 
back and processing leaves for the fungus; there are (1 — q)pA workers taking care of the fungus 
garden and there is healthy fungus F in the garden. Thus, fungus growth is a product of two 
different sets of tasks performed by workers which can be represented by the following diagram: 



qpA + {l — q)pA +F —> F + new fungus 

energy from workers collecting and processing leaves energy from workers tending fungus 

Therefore, by applying the concept of the kinetics of functional response (Real 1977) we can assume 
that the numerical response of fungus to ants is a Holling Type III function 

b + p^q{l- q)A^ 

where h is the half-saturation constant. The population of fungus decreases due to the consumption 
by ants and its mortality. Here, we assume that the fungus suffers from density-dependent mortality 
due to self-limiting (Holland and DeAngelis 2010). Thus, the population dynamics of fungus can 
be described by follows: 

dF _ (r,p\il-q)A^ -d,F-racA\F (2) 



dt \b + p'^q(l-q)A'^ 

where is the maximum growth rate of the fungus; c is the conversion rate between fungus and 
ants and d/ is the mortality rate of the fungus. See Table [l] for the biological meanings of the 
completed list of parameters. 

Remark: Holling (1959) and Murdoch (1969) have discussed the application of a type HI func- 
tional response associated with learning processes of predators being able to adjust their feeding 
rate actively based on the quantity and density of available prey. An ant colony similarly changes 
allocation of workers to different tasks as needs change with changing colony size (HoUdobler and 
Wilson 1990). In this way the ants actively modify their time expenditure and task allocation just 
as predators actively modify their feeding behaviors. 
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Table 1: Biological meanings of parameters in the system l[T]l-ll2^ 



Parameters 


Biological Meaning 


ra 


Maximum growth rate of ants 




Maximum growth rate of fungus 


c 


Conversion rate between fungus and ants 


da 


Death rate of ants 


df 


Death rate of fungus 


P 


Proportion of ants that are workers 


q 


Proportion of workers that take care of fungus 


b 


Half-saturation constant 



3. Mathematical Analysis 

Let a = p^q{l — q) and Tc = cva, then based on the assumptions listed in Section 2, an interaction 
between ants and fungus at the early stage of colony may be modeled by the following differential 
equation: 

(7/1 

^ = {raF-daA)A (3) 

where a can be considered as a parameter measuring the division of labor in the colony of ants and other 
parameters are strictly positive. Since 

a = p^q[l-q), p€ [0, 1] and q € [0, 1], 

therefore, 

a e [0, 0.25] 

where a achieves its maximum 0.25 when p ~ ^,q — 1/2. In reality, p is always less than 1 since the 
biomass of queen and other stages of ants (i.e., larvae, pupae) other than adult ants is greater than 0. 
Thus, a is strictly less than 0.25 in the real biological system. 

Lemma 3.1. The system ([Sjl-Q is positively invariant and hounded in M^. In particular, if both 
A[Q) > and F{0) > 0, then A{t) > and F{t) > for all t > 0. 

Proof Notice that A{0) = then A{t) = for all t > 0; F(0) = then F{t) = for aU t > 0. Therefore, 

1. If yl(0) = and F(0) = 0, then {A{t), F{t)) = (0, 0) for all t > 0. 

2. If A{0) = and F(0) > 0, then 

dF ( rtaA^ , A , , 

^ - dfF -rcA] F = -dfF^ < 0. 



dt \b + aA^ 
Thus, limt^oo F{t) = hmt^^o = 0. 

3. If F(0) = and A(0) > 0, then 

dA 



{raF ~ daA) A = -daA' < 0. 
dt 



Thus, limt^oo A{t) = limf_^oo = 
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If A{Q) > and F{0) > 0, then due to the continuity of the system, it is impossible for either A{t) or 
F{t) to drop below 0. Thus, for any A(0) > 0, F(0) > 0, we have A{t) > and F{t) > for all t > 0. 
Now assume A{0) > 0,F(0) > 0, then according to the expression of ^^Jf^ , we have 



dF 



TfaA^ 
b + aA'^ 



dfF -rcA] F < {rf - dfF) F. 



Thus, limsup(_j.oo F{t) < 



gy. This indicates that for any e > 0, there exists T large enough, such that 



e for all t > T. 



Therefore, we have 



^ = (r,F - daA) A < (^aCj-^ + e) - daA ] A, for aU t > T. 



Since e can be arbitrarily small, thus limsup^^g^ A{t) < -^f^- Therefore, we have shown that the system 
(|3|-([4| is positively invariant and bounded in R^. More specifically, the compact set [0, ^f^] x [0, ^] 
attracts all points in Mi. 



Moreover, if both A{0) > and F(0) > 0, then we have follows 



dA 
lit 
dF 

~dt 



{raF ~ daA) A > -daA^ 

rjaA^ 



A(t)> 



m 

l + dj 



> 



b + aA'^ 



dfF -TcA] F> -dfF^ ^ F{t) > 



m 

1 + dft 



> 



Therefore, if both A{0) > and F{0) > 0, then A{t) > and F{t) > for aU t > 0. 



□ 



Remark: Lemma |3 . 1 1 indicates that the population of leaf cutter ants and their fungus is boun ded due 



3.1 



to the limited resource in nature. Let = {{A,F) £ ; yl > 0, > 0}, then from Lemma 
know that (A(0), F(0)) G implies that {A{t), F{t)) e M.^ for aU t > 0. 

Proposition 3.1. For any initial condition taken in M'^, the trajectory of the system (|3])-(|4| is converg- 
ing to an equilibrium point. 



Proof. By Poincare-Bendixson Theorem (Guckenheimer and Holmes 1983), the omega limit set of the 
system Q-Q is cither a fixed point or a limit cycle. If there exists a function B{A,F) : M"^ — > R 
such that 



±[B{A,F){raF-daA)A] + ^ 



B{A,F) 



rfoA^ 
b + aA^ 



dfF-r,A\F 



< 0, 



then we can use Dulac's criterion (Guckenheimer and Holmes 1983) to exclude the existence of a limit 

1 

AF- 



cycle for the system . Let B{A,F) = Then, 



d_ 
dA 



[B{A, F) {raF - daA) A] 



d_ 
OF 



B{A,F) 



rjaA^ 
b + aA^ 



-dfF 



A F 



da 

' F 



A 



< 



holds for any {A{0),F{0)) £ M^. Therefore, by Dulac's criterion, the system ([sjl-Q has no limit cycle, 
i.e., any trajectory of ([3|-([4]) starting with a non-negative initial condition converges to a fixed point. □ 



Remark: Proposition |3.1| implies that the population dynamics of leaf cutter ants and their fungus is 
simple in the sense that they do not have a limit cycle, i.e., if time t is large enough, then the population 
of leaf cutter ants and their fungus approach to some fixed point. Therefore, the short time dynamics 
of leaf cutter ants and their fungus garden is more important since it can give us more information on 
the dynamics of the interaction between ants and fungus. 
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Proposition 3.2. If a < Ah ^ , then the system has only trivial equilibrium (0,0); 

while if a = ''" '^^J ^ , then the system (|3|-(|4]) has the only positive equilibria 

rj Ta rf ra da 



{A\F') = 



2 {rc ra + df da) ' 2ra {rc ra + df da) 



in addition to (0, 0); while if a > Ab ^ ^° ^ , then the system ([Sj-Q /las the following two positive 
equilibria in addition to (0,0); 



{A'\F'') = (^A'\ ^A'l^ and {A'^,F'^) = (^A^^, ^A^^^ 



(5) 



where 

A' 



Tfra ( rfra \^ ^ 



'2{rcra + df da) y V 2 (f c ra + df da) J a 



^ "Tfra , J( rfra 



2{rcra + dfda) \l \2 {r^ra + df da) 



2 



b 



Proof. It is easy to see that (0, 0) is always an equilibrium of the system Q-Q. The nullclines of (|3])-(|4| 
can be founded as 

dA , da , 

— = 0^A = or F ^ —A 
dt ra 

dF rtaA^ r^ , 

-0^F = or F = —I -r- - -^A 



dt dfib + aA"^) df 

df{b+aA^) ^ Tf ' 



By solving , a2\ — = for A, we have the following two cases 



1. If a > Ab i ^" '"^^^'^ ^ , then by simple algebraic calculations, there are the following two positive 
solutions of a^^g^ ~ 37^ = 



A' 



rfra j ( rfra \^ ^ 



2{rcra + dfda) V \2 {r^ra + df da) J a 



^ rfra , J( rfra 



2{rcra + df da) ^ \2 (r^ ra + df da) 
Thus, the two interior equilibria are 



2 



b 



{A'\F'') = [A'\^A'^j and {A'^,F'^) = f A'^^A^^ 



2. If a = 46 ^ ^^1^ j , then the system (|3])-(|4| has only one positive equilibria F*) where 

rfra rfra da 



{A\F' 



2 {rc ra + df da) ' 2ra {rc ra + df da) 
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3. If a < 46 ( j"^!^'^^^ "^" ^ , then there is only one trivial equilibrium: A — and F = 0. 



Therefore, the statement of Proposition |3.2| holds 



□ 



Remark: Recall that a is a parameter measuring the division of labor of workers [j 9]. Proposition 3.2 
implies that if a is too small, i.e., the ratio of adult ants that take care of fungus to adult ants that 
forage for leaves, is too small, then the system ([3|-(|4]) has only trivial equilibrium point (0, 0). This 
leads to the following theorem: 

Theorem 3.1. [Extinction of Two Species] If a < 4b ( j' ^^^^f ) ' then the system has 
global stability at (0,0). 



Proof. From Proposition 



3.2 



we know that for any initial condition taken in M^, the trajectory of 



the system ([3|-([4]) is converging to an equilibrium point. If a < 45 ^^^^^^r-^^Jj-^^ , then according to 
Proposition |3.2[ the only equilibrium of the system ([3|-(|4| is the origin (0, 0). Therefore, we can conclude 

that the system (|3])-(|4| has global stability at (0, 0) when a < 45 ^ — 



ra-\-df da 



□ 



Biological Implications: Theorem 3.1 indicates that division of labor is an important factor deter- 
mining whether the early colony stage of leaf cutter ants can survive or not. Recall that the proportion 
of ants performing a task is essentially equivalent to energy devoted to a given task. In the case that the 
population of adult ants is too small, i.e., q is too small, or the population of adult ants foraging for food 

is too small, i.e., (1 — q) is too small, then a < q{l — q) will be too small such that a < 45 ( j" ^!]^^^ ^ . 
This leads to the extinction of both ants and fungus. 

In order to investigate the biological conditions when leaf cutter ants and their fungus can coexist, 
we have the following theorem: 

Theorem 3.2. [Coexistence of Two Species] If a > 45 ( j" ^'[^'^1 ^° ^ , then the system ([sjl-Q has 

two positive equilibria and where is always unstable and (A'^,^'^) is 

always locally asymptotically stable. 



Proof. From Proposition 3.2 we know that the system ([Sj-Q has two positive equilibria (A*^, i^*^) and 

(^'^, F*^) when a > 45 { ^^ ^!]^^J ^ . The local stability can be determined from the eigenvalues of its 
Jacobian Matrices evaluated at these equilibria. 

Assume that {A* , F*) is an equilibrium point of ([sjl-Q, then its Jacobian Matrices evaluated at this 
equilibrium can be written as follows 



\{A',F') 



dr. df 



-da A* 
I- A* 



TaA* 



(6) 



Then we have 



trace {J\(A',F')) 
det {J\(A*,F')) 



-dr, A 



Tg+df 



< 



da {A*Y{ra ra+da df){2ra Ta A'-Tg rf+2da df A*) 
df 



This implies that if A* > 



2{ra r^+da df) 



then {A*,F*) is locally asymptotically stable; while if A* < 
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2(r r'^+d dj) ' then {A* , F*) is a saddle node, i.e., unstable. Since 



"^{rcra + df da) \l \2 {rcVa + df da) J a 2{rarc + dadf) 




A- = + J( "^l^ ] --> '"'f 



2{rcra + dfda) V \2 (rcTa + df da) J a 2{rarc + dadf) 



Therefore, is always unstable and is always locally asymptotically stable when 

Biological Implications: Theorem |3.2| implies that if allocation of workers to different tasks is in a 
good range, i.e., a > 46 ( j" '"'^^^^ ^ , then both leaf cutter ants and their fungus garden can coexist, be- 
cause the system (|3|-([4| has a locally asymptotically stable interior equilibrium (A*^, F*^) . On the other 

hand, for a fixed value of a, if da, rc and are small enough, then a > 4& ( j" ^"^^^ ^ holds, thus two 
species can coexist. Now the more interesting question is whether relative allocation among tasks is the 
only factor determining whether ants and fungus can coexist. The next theorem will answer this question. 



Theorem 3.3. [Basin of Attraction of (0,0)] The trivial equilibrium (0,0) is always locally asymp- 
totically stable if a ^ 4b (^^-^-^-f-^-^j-^^ ■ Moreover, if a > Ab ^!]i~^^^ ^ , then the basin of attraction 
of (0, 0) contains in the region i?(o,o) \ 'S'(a»i._f»i) where 



and 



S^A^.^pn. = {{A,F) e Rl : hm {A{t),F{t)) = {A^\F^^)}. 



Proof. If a < 46 ( j' ^r^rf ) ' then according to Theorem js.lj (0,0) is global stable in M^, thus it is 



locally asymptotically stable. Now we need to consider the case that a > 46 
First, we claim that the region defined by 

f . . X "9 rtaA^ rrA daA 

rti = \ {A, F)eRl: _ ^ < ^ < ^ 

\ + df[b + aA^) df - - ra 

is positively invariant. Assume that this is not true. Then there is some initial condition (A(0),-F(0)) 
take n in fii such that for some future time T such that {A{T),F{T)) is leaving fli. From Proposition 
we know that the system ([sjl-Q has only one equilibrium point in Qi, thus due to the 



3.2 



continuity of the system, there exists some time T such that we have one of the following two cases: 
1. For aU < i < T, 

rfoA^t) r,A{t) ^ ^^^^ ^ daA{T) ^^ 



dfib + aA^t)) df 

at t = T, 



r^aA^iT) ^c^(^) _ ^(^) F(r) < 



dfib + aA^T)) df 
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and for some e > and T < t < T + e, we have 

rfaA'^it) rcA{t) daA{T) 



> F{t) and F{t) < 



2. For all < t < T, 



at t = T, 



df{b + aA^t)) df 

df{b + aA^t)) df ^ Ta ' 

TfaA^iT) r,A{T) ^ , ^^rj.. _ daA{T) 



df{h + aA^{T)) df 
and for some e > and T < i < T + e, we have 

daA{t) 



< F{T) and F{T) = 



nt)> ^ 

' a 



If the first case holds, then at time t = T we have 
dA, 



^\t=T = iraF{T)-daA{T))A{T) = 

2 



dF _ ( rfaA^{T) 



This implies that there exists some small e such that A{t) < A{T), F{t) < F{T) for all T < t < T + e, 
which contradicts the conditions for the first case. Similarly, we can show it is impossible for the second 
case to be held. Therefore, fli is positively invariant. 

Now we will show that -B(o,o) is positively invariant. Define 

1^2 = -8(0,0) \ ^1 and = f^i \ {{A'\F'^)} . 

Then, fis is also positively invariant since (A*^, F*^) is an equilibrium point and -B(o,o) = U ^3. For 
any initial condition {A{0),F{0)) taken in i?(o.o)i there are the following two cases: 

1. If (^(0),F(0)) e n-ii, then {A{t),F{t)) e O3 for alH > since O3 is positively invariant; 

2. If (A(0),F(0)) e ^2, then either {A{t),F{t)) G ^2 for alH > or there exists some T such that 

^\t=T = {raF{T)-daA{T))A{T) = 

dF _ f rfaA^iT) 



dt 1*=^ 



This implies that {A{T),F{T)) e O3. Since O3 is positively invariant, then {A{t),F{t)) e fls for 
all t > T. 

Therefore, -B(o,o) is positively invariant. 

Define a Lyapunov function V = A"-F^ : -B(o,o) ~^ where both a and (3 are positive. Then we 
have 



= c^A'^-'F^^ + pA^F^-'^^ = aA"F^ {r^F - d^A) + pA-^F^ {j^^ ' ' 



= V 
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Choose a,/3 such that — — 0, i.e., ^ = gy, then the expression of ^ can be simphfied as 



dV 
Itt 



= V 



dfda 



A 



= 13V A 



rfaA-b[r, + ^^j -a{r, + ^]A^ 
b + aA"^ 



Define f{A) = rjaA — b + ^^^j — a yvc + ^^^^r^j A^, then the sign of ^ depends on the sign of 

f{A). Notice that f{A) = {A- A''^){A'-'^ - A) is negative if < A < A'^ and any point {A, F) e S(o,o) 
satisfying < ^ < < F < F*i. 

Since -B(o.o) is positively invariant, for any initial condition taken in _B(-o_o)j have ^ < for 
all future time. This indicates that A{t) and F{t) approach to some fixed point contained in 5(0,0) 
(the closure of i?(o,o))- Notice that -B(o,o) contains only (0,0) and (v4*^,F'^). According to Theorem 



3.2 



(^*^,F*^) is unstable when a > 4b ( j'^^^^J '^" '^ ■ Then based on Hartman-Grobman Theorem 

(Robinson 1998), any point in i?(o.o) \ 'S'(a»i,f*i) will not approach to (A*^,F*^), and therefore, it will 
approach to (0, 0). 

Therefore, the statement of Theorem 13.31 holds . □ 



A direct corollary of Proposition |3.1[ Theorem |3.2| and Theorem |3.3| is as follows 



Corollary 3.1. If a > 4b 



Ta+df da 



then the system has two attractors 

(0,0) and {A''^,F'^) . 



If the initial condition (^(0), F{0)) is too small such that it contained in 5(o,o) \ ^(A' 

hm (A(t),F(i)) = (0,0); 

while the initial condition {A{Q)^F{Q)) is large enough, then 

hm {A{t),F{t))^{A^\F^'). 



the 



Biological Implications: Theorem |3.3| and Corollary |3.1| suggests that the initial population of leaf 
cutter ants and fungus is another important factor that determines whether ants and fungus can coexist 
or not. If initial population is contained in -B(o,o) \ •S'f^ii.i^ii) , then both ants and fungus will go extinct 

even if the division of labor is in a good range, i.e., a > 4b (^ ^'' ^ . 



4. Numerical simulations, data and sensitive analysis 

In this section, we validate our model Q-Q by performing numerical simulations, sensitivity analysis 
and parameter estimations based on the experiment data. The numerical simulations fit the data very 
well (see Figure [l] and [2]) , which suggests that our model ([Sjl-Q is well defined. Sensitivity analysis 
around these chosen parameter values provides information on the governing factors for the ecological 
process modeled by ([sjl-Q. Parameter estimations with different initial guessing values indicate that 
the population dynamics of leaf cutter ants and their fungus may be unstable at the early stage of colony 
expansion, as supported by the empirical data (Clark and Fewell in preparation). 



11 



Numerical simulations and experimental data 

In this subsection, we compare the numerical simulations of the model (|3|-(|4]) by using parameter 
values in certain intervals (see these values in Table [2]). These intervals are obtained from the approxi- 
mations according to data and literature (Brown, Bot and Hart 2006; Clark and Fewell in preparation). 
Table[2]lists the range of parameters and the specific values (i.e., ra = 0.1, r/ = .7, da = 0.1, d/ = 0.2,6 = 
0.002, rc = .0045,0 = 0.2,^(6) = 0.05, F(6) = 0.3) that generate dashed lines in Figure [l] and [2] from 
week 6 to week 29. Other values in the interval can generate similar dynamics as the chosen values, but 
the chosen values fit the data well. 



Table 2: Intervals and chosen values of parameters in the system l|3]l-H 



Parameters 


Intervals 


Chosen values 


ra'. Maximum growth rate of ants 


(0.05, 0.3) 


0.1 


r f : Maximum growth rate of fungus 


(0.01,1) 


0.7 


rc'. Conversion rate between fungus and ants 


(0.001,10) 


0.0045 


da- Death rate of ants 


(0.001,1) 


0.1 


df] Death rate of fungus 


(0.001,1) 


0.2 


b: Half-saturation constant 


(0.001,10) 


0.002 


a: Measurement of the division of labor 


(0, 0.25) 


0.2 


yl(6):Biomass of ants at week 6 


(0.001, 0.1) 


0.05 


i^(6):Biomass of ants at week 6 


(0.001, 1) 


0.3 



Figure [T] and Figure [2] provide the comparison between ecological data (solid lines with error bars) 
and simulations (dashed lines) generated by the model ([Sjl-Q when = 0.1, = -7, da = 0.1, rf/ = 
0.2, b = 0.002, rc = .0045, a = 0.2, A{6) = 0.05, F(6) = 0.3. 

1. In Figure [T] the left figure is the biomass of Ants v.s. Time in weeks and the right figure is the 
biomass of Fungus v.s. Time in weeks. By comparison, we can see that the simulations fit the data 
very well, especially for biomass of fungus. Overall, the simulation of the biomass of ants is larger 
than the experimental data (the right figure). This is expected, because the equation ^ models 
the biomass of all ants including the queen, eggs, larvae, pupae and workers, while the experiment 
only measures the biomass of workers. In addition, both data and simulations suggest that ants 
have exponential growth while the fungus has linear-like growth from weeks 6-29. Recall that our 
focus is the ergonomic growth stage of the ants which starts when the first workers appear. The 
exponential growth of ants at this growth stage confirms the study by Oster and Wilson (1978). 

2. The right figure of Figure[2]represents logiQ(the biomass of ants-|-l) v.s. logiQ(the biomass of fungus-|-l), 
which provides the information on the relationship between the growth rate of ants and the growth 
rate of fungus. Simulations fit data extremely well. Both suggest that the growth rate of ants and 
fungus increase over time, and the growth rate of ants increases faster than fungus, which may be 
caused by changes in the efficiency of the conversion between ants and fungus at the early colony 
stage. 

The left figure of Figure[2]is the ratio of log]^Q(the biomass of ants-t-1) to log]^o(the biomass of fungus-f-1) 
v.s. time in weeks, which provides information on the relative growth rate of ants to fungus: the 
simulation fits data very well from week 10 to week 29 but shows some inconsistency between the 
data and the model fitting during week 6 to week 9. In this case the model is a more accurate 
descriptor of population dynamics than the collected data because the biomass of the ants during 
this time consists almost entirely of immature workers or ants in the larvae/pupae stage. The 
data do not account for this ant biomass and thus, from week 6 to week 9 the ant population 
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Figure 1: The solid lines with error bars are from data and the dashed lines are simulations generated from the model 
when ra = 0.1, = .7, da = 0.1, dj- = 0.2,6 = 0.002, Tc = .0045, a = 0.2,71(6) = 0.05, F(6) = 0.3. The left figure is 
the biomass of Ants v.s Time in weeks and the right figure is the biomass of Fungus v.s. Time in weeks. 



may be largely underrepresented. This under-representation leads to an increase in the slope of 
data when the actual result should be closer to the model output during this time. Thus, the 
possible explanations for the inconsistency between the data and the model fitting during week 6 
to week 9 can be summarized as follows: 1. The equation ^ models the biomass of all ants while 
the experiment only measures the biomass of workers; thus, the model ([sjl-Q should generate the 
larger ratio of log^oCthe biomass of ants+1) to log]^Q(the biomass of fungus+1); 2. For the first few 
weeks (week 6-9), the real population dynamics of ants and fungus are highly unstable and may 
have very different ecological properties than our model assumptions. Stochastity and multiple life 
stages of leaf cutter ants (e.g., eggs, larvae, pupae) may be considered in future models. 

The fitting of the model to the data for the ants and fungus population is evident from Figures [l] 
and [2] The comparison between simulations generated by the model and data suggests not only 

that the fit of the model to data is accurate but also that parameters match with expected values for 
growth, death, and division of labor. A recent experiment study by Clark and Fewell (in preparation) 
on leaf-cutter ants shows that the parameter values for ra,rf,df that generate Figures [l] and [2] are very 
close to actual data. The initial condition used for generating Figures [l] and [2] is the mean value from 
the experimental data by Clark et al (preprint). A study by Brown, Bot and Hart (2006) on mortality 
rates of leaf cutter ants and division of labor suggests that the death rate matches with the parameter 
da- Notice that the values of parameters such as rc,b are difficult to measure in the experiments. The 
good fit of the model to data (see Figures [T] and [2]) generated by the values listed in Table [2] provides an 
approximation of and b. In the following subsection, we will examine the sensitivity of these parameter 
values and the initial condition. The parameter d/ is difficult to measure experimentally because of the 
efficiency of the mutualistic relationship between the ants and fungus. 

The leaf-cutter ants cultivate and consume the fungus inside underground gardens, so data collection 
in the field ultimately results in the death of the colony. Our data measures of fungal biomass include 
estimates of fungal loss from waste piles; however, calculating plant waste underestimates the death rate 
of the fungus, because much of the fungus is consumed or used for other purposes inside the colony. By 
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Figure 2: The solid lines with error bars are from data and the dashed lines are simulations generated from 
the model when = 0.1, rj = .7, da = 0.1, d/ = 0.2,6 = 0.002, = .0045, a = 0.2, yl(6) = 

0.05, F{6) = 0.3. In the left figure, the x-axis represents time in weeks and the y-axis represents the ratio of 
logj^Q(thc biomass of ants+1) to logjg(the biomass of fungus+1); while in the right figure, the x-axis represents represents 
the value of logj^Q(thc biomass of ants+1) and the y-axis represents the value of logjQ(the biomass of fungus+1) 



the time the fungus biomass has been disposed of, it is only a fraction of its original biomass, making 
any measurement of its mass/area inaccurate. AdditionaUy, actual consumption of the fungus cannot 
be easily observed in a field or laboratory setting. 

4-2. Sensitivity analysis 

Input factors for our mathematical model ([sjl-Q consist of seven parameters and two initial con- 
ditions for independent and dependent variables of the model. Because of natural variation, error in 
measurements, or simply a lack of current techniques to measure some parameters, it is necessary to 
perform sensitivity analysis to identify critical inputs (parameters and initial conditions) of our model 
and quantifying how input uncertainty impacts model outcomes (i.e., the dynamics of the ants and fun- 
gus biomass of ants A{t),F{t)). In this subsection, sensitivity measure of the model (|3])-(|4| is computed 
numerically by performing multiple simulations varying input factors around the nominal value listed in 
Table [2 

Sensitivity analysis for parameters and initial conditions were performed using an extension of the 
MATLAB function ODE23tb, a stiff solver for ordinary differential equations. The extension main- 
tains the same calling sequences as ODE23tb; it differs only in the algorithm used to compute the 
derivatives. The algorithms used in this instance are the internal numerical differentiation and it- 
erative approximation based on directional derivatives methods described by H.G. Bock (1981) and 
T. Maly with L. R. Petzold (1996) respectively. The output of the function is similar to that of 
ODE23tb with an additional array containing the derivatives (sensitivities) of the solution with re- 
spect to a given parameter vector. The sensitivity of all parameters and initial conditions around 
r„ = 0.1, r/ = .7, da = 0.1, df = 0.2,6 = 0.002, = .0045, a = 0.2, A(6) = 0.05, i^(6) = 0.3 from week 6 
to week 29 are shown in Figure |3][7) 
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Figure 3: The left figure is the sensitivity of the maximum growth rate of ants Ta and the right figure is the sensitivity 
of the maximum growth rate of fungus rf for the model l|3]l-(|4]l around the chosen values where Va = 0.1, = .7, da = 
O.l, df = 0.2, b = 0.002, = .0045, a = 0.2, A{6) = 0.05, F(6) = 0.3. 



As we would expect, the effects of all the parameters are strictly cumulative. The effects of each 
parameter at the beginning of the colony are relatively small. As time progresses the parameters have 
a much greater effect on the model, most of them having their largest effect on the biomass of fungus 
at around week 25. There is a time shift between the effect on the biomass of ants and the biomass of 
fungus: the largest effect on the biomass of fungus is always earlier than the largest effect on the biomass 
of ant. In addition, all the parameters have larger effects on the biomass of fungus than on the biomass 
of ants. The detailed results on sensitivity analysis for each parameter and the initial condition can be 
summarized as follows: 

1. Figure [3] shows that the growth parameters and rf have a positive effect on the model as a 
function of time; as the model progresses, and rj have a larger effect on the increase of both 
biomass of ants and fungus. The largest effect of both and on biomass of fungus occurs at 
week 25, while the effect of both Va and on biomass of ants is an increasing function of time 
and achieves the largest effect at the end of experiment, i.e., week 29. 

2. Figure |4] shows that the death parameters da and df have a negative effect on the model as a 
function of time. As the model progresses forward in time the decrease in the biomass of ants and 
fungus respectively is increasingly affected by parameters da and df. The largest effect of and 
r/ on biomass of fungus occurs at week 25 and week 26 respectively while the effect of both da 
and df on biomass of ants is an increasing function of time and achieves the largest effect at the 
end of experiment, i.e., week 29. 

3. Figure [6] shows the conversion rate between ants and fungus rc has a negative effect on the output 
of the model while the measurement of the division of labor of ants has a positive effect. The 
effect of Tc on the biomass of both ants and fungus is a decreasing function of time while the 
division of labor parameter a shows very similar behavior to the growth parameters and ry; 
this suggests that by maximizing the efficiency of division of labor, both the ants' and fungus's 
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Figure 4: The left figure is the sensitivity of the death rate of ants da and the right figure is the sensitivity of the death 
rate of fungus df for the model l|3]l-(|4]l around the chosen values where = 0.1, = .7, da = 0.1, df = 0.2, b = 0.002, Tc = 
.0045, a = 0.2, A{6) = 0.05, F(6) = 0.3. 



population growth will be maximized. This agrees with article main purpose number 2. 

4. Figure [t] shows that both initial condition A{6) and F{6) have positive effects on the model as a 
function of time. The effect of A{6) on biomass of both ants and fungus is an increasing function 
of time. The effect of F{6) on the biomass of ants is an increasing function of time while the effect 
on the biomass of fungus is decreasing until week 15 and then increasing until week 29. 

5. Figure [5] shows that the half-saturation constant b has a negative effect on the model as a function 
of time. As the model progresses forward in time, the decrease in biomass of ants and fungus 
respectively is increasingly affected by b. The largest effect of 6 on fungal biomass occurs at week 
25 while the effect of b on biomass of ants is an increasing function of time and achieves the largest 
effect at the end of experiment, i.e., week 29. Notice that b has the largest sensitivity among all 
the parameters and initial conditions. 

In conclusion, the growth parameters ra,rf, the division of labor parameter a and the initial con- 
ditions A{6),F{6) have positive effects on the biomass of ants and fungus while the death parameters 
da,df, the conversion rate and the half-saturation constant b have negative effects on the biomass of 
ants and fungus. Among all these parameters and the initial conditions, the parameter b has the largest 
sensitivity and the conversion rate rc, the death rate of ants da and the initial value of fungus F{6) have 
relative small sensitivity to the output of the model. 

4.3. Parameter estimations 

In this subsection, we use Nonlinear grey-box models (System Identification Toolbox provided in 
MATLAB) to perform parameter estimations based on experimental data by fixing the estimated interval 
of a to be (0, 0.25) and the estimated intervals of other parameters to be (0, 00) . The estimated values can 
be varied depending on the initial guesses and the estimated intervals. Table [3] lists the best estimated 
values and their standard deviations when the initial guessing are the values that generated Figures [T] 
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Figure 5: The sensitivity of the half-saturation constant b for the model l|3]l-l[4]| around the chosen values where = 
0.1, = .7, da = 0.1, df = 0.2,b = 0.002, = .0045, a = 0.2, A(6) = 0.05, = 0.3. 



and [2j Table |4] lists the best estimated values and their standard deviations when the initial guesses are 
different from the one 's listed in Tabled 



Table 3: Parameter estimations of the system JSjl-H case one 



Parameters 


Initial values 


Intervals 


Estimated values 


Standard Deviation 




0.1 


(0,oo) 


0.0436 


0.0949 




0.7 


(0,oo) 


1.1667 


4.3344 




0.0045 


(0,oo) 


4.9781e-12 


3.5671 


da 


0.1 


(0,oo) 


1.5144e-07 


0.6700 


df 


0.2 


(0,oo) 


0.2707 


0.3035 


b 


0.002 


(0,(X3) 


0.0050 


34.2259 


a 


0.2 


(0,0.25) 


0.1656 


1140.99 


A{6) 


0.05 


(0,c^) 


0.0988 


DNE 


F(6) 


0.3 


(0,oo) 


0.1860 


DNE 



The comparison between Table |3] and Table |4] can be summarized as follows: 

1. The different initial guesses values will give different estimated values. This may be caused by the 
fact that the model ([3|-(|4| has multiple attractors. 

2. The estimated values of rc and da are both extremely small. In addition, the smaller the initial 
guesses of da and r^, the smaller the estimated values of these parameters. This may suggest that 
Vc and da have little effect on the population dynamics of ants and fungus at the early colony stage, 
which has been confirmed by their small sensitivity (see Figure H and pi) . 
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Figure 6: The left figure is the sensitivity of the conversion rate between ants and fungus Tc and the right figure is 
the sensitivity of the measure of the division of labor of ants a for the model |3]l-l[4]| around the chosen values where 
Ta = 0.1, = .7, da = 0.1, d/ = 0.2,6 = 0.002, = .0045, a = 0.2, A(6) = 0.05,F(6) = 0.3. 



3. The standard deviations of estimated a, 6, r^, da are extremely large, which may be caused by two 
factors: 1. These parameters are not independent; 2. The extremely small value of da and r^. 

The summary above indicates that the population dynamics of ants and fungus may be highly unstable 
at the early stage of colony development. Notice that collected data is from the successful colonies only. 
The extremely small estimated value of and da may suggest that the conversion rate between ants 
and fungus and the death rate of ants are not as important as other factors such as the growth rate 
parameters ra,r/ and the death rate of fungus df. Possibly, a multiple-stage model that includes the 
stages of eggs, larvae, pupae or even a stochastic model should be introduced in order to get a better 
understanding of the detailed ecological processes. 

5. Conclusion 

In this article, we develop a simple mathematical model (|3|-(|4| to study mutualism interactions 
between leaf cutter ants and their fungus garden at the early colony stage with the following unique 
features 

1. The net benefit of the obligate fungus to leaf cutter ants is determined by the difference between 
the overall performance of collecting leaves and cultivating fungus by worker ants and the amount 
of fungus eaten by queen, larvae and workers; while the net benefit of obligate ants to fungus is 
determined by the difference between the amount of consumed fungus and the mortality rate due 
to the energy spent on collecting leaves and cultivating fungus. 

2. The division of labor of leaf cutter ants: workers perform different tasks to maintain their fungus 
gardens. This feature allows us to apply the concept of the kinetics of functional response to model 
the numerical functional response of fungus. 
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Figure 7: The left figure is the sensitivity of the biomass of ants in week 6 A{6) and the right figure is the sensitivity of the 
maximum growth rate of fungus Vf for the model l|3]l-(|4]l around the chosen values where = 0.1, ry = .7, da = 0.1, df = 
0.2, b = 0.002, rc = .0045, a = 0.2, A{6) = 0.05, F(6) = 0.3. 



The mathematical analysis of ([Sj-Q gives the completed global dynamics of the model (Theorem 
3.1 3.2 3.3 and Corollary 3.1). These theoretical results suggests that: 1. The division of labor of ants 



can determine whether leaf cutter ants and their fungus garden are able to coexist; 2. When the division 
of labor is in a good range, the initial populations of leaf cutter ants and fungus are needed to be larger 
than some threshold in order to coexist. 

Finally, we validated the model Q-Q using empirical data. The comparison between model sim- 
ulations and data supports the fact that ([Sjl-Q is well defined for modeling the population dynamics 
of the leaf cutter ants and fungus during the incipient colony stage (the early ergonomic growth stage). 
The good fit between the model and data also provides us an approximation of the values of difficult 
measured parameters such as the conversion rate between fungus and ants and the half-saturation 
constant b. Sensitivity analysis implies b has the largest effect on the output of the model. Both sensi- 
tivity analysis and parameter estimations suggest that the population dynamics of early stage colonies 
may not be stable and the growth rate parameters ra,rf and the death rate of fmigus df are important 
factors for determining the population dynamics for the successful colony. 

The inconsistency between the data and the model fitting during week 6 to week 9 (Figure [2| suggests 
that a more realistic and detailed model is needed during this period. Thus, consideration of multiple life 
cycle stages, including eggs, larvae and pupae, or even stochastity should be included in further modeling 
work. In addition, for future experiments, the biomass of, larvae and pupae should be measured as well. 
These different life cycle stages are likely to represent a larger proportion of the total ant biomass in 
early stages of colony growth and additionally likely have more variance. This could be our future work. 
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Table 4: Parameter estimations of the system Jsl-H case two 



Parameters 


Initial values 


Intervals 


Estimated values 


Standard Deviation 




0.075 


(0,oo) 


0.0747 


0.0230 




0.15 


(0,oo) 


0.1585 


0.2055 




0.0001 


(0,oo) 


9.9447e-05 


0.3906 


da 


0.0001 


(0,oo) 


3.1570e-07 


0.5807 


df 


0.03 


(0,oo) 


0.0292 


0.1061 


b 


0.00001 


(0,oo) 


8.2187e-05 


33.6161 


a 


0.2 


(0,0.25) 


0.2004 


81819.3 


A{6) 


0.05 


(0,oo) 


0.0498 


DNE 


m 


0.3 


(0,oo) 


0.2995 


DNE 
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